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LISA Pathfinder (LPF) is a mission aiming to test the critical technology for the forthcoming 
space-based gravitational wave detectors. The main scientific objective of the LPF mission is to 
demonstrate test-masses free-falling with residual accelerations below 3 x 10~^* m s^'^/v'Hz at 1 
mHz. Reaching such an ambitious target will require a significant amount of system optimisation 
and characterisation, which will in turn require accurate and quantitative noise analysis procedures. 
In this paper we discuss two main problems associated with the analysis of the data from LPF: 
i) Excess noise detection and ii) Noise parameter identification. The mission is focused on the 
low frequency region ([0.1, 10] mHz) of the available signal spectrum. In such a region the signal 
is dominated by the force noise acting on test masses. At the same time, the mission duration 
is limited to 90 days and typical data segments will be 24 hours in length. Considering those 
constraints, noise analysis is expected to deal with a limited amount of non-Gaussian data, since 
the spectrum statistics will be far from Gaussian and the lowest available frequency is limited by 
the data length. In this paper we analyze the details of the expected statistics for spectral data 
and develop two suitable excess noise estimators. One is based on the statistical properties of 
the integrated spectrum, the other is based on Kolmogorov-Smirnov test. The sensitivity of the 
estimators is discussed theoretically for independent data, then the algorithms are tested on LPF 
synthetic data. The test on realistic LPF data allows the effect of spectral data correlations on the 
efficiency of the different noise excess estimators to be highlighted. It also reveals the versatility of 
the Kolmogorov-Smirnov approach, which can be adapted to provide reasonable results on correlated 
data from a modified version of the standard equations for the inversion of the test statistic. Closely 
related to excess noise detection, the problem of noise parameter identification in non-Gaussian data 
is approached in two ways. One procedure is based on maximum likelihood estimator and another 
is based on the Kolmogorov-Smirnov goodness of fit estimator. Both approaches provide unbiased 
and accurate results for noise parameter estimation and demonstrate superior performance with 
respect to standard weighted least-squares and Ruber's norm. We also discuss the advantages of 
the Kolmogorov-Smirnov formalism for the estimation of confidence intervals of parameter values 
in correlated data. 
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I. INTRODUCTION 



LISA Pathfinder (LPF), a European Space Agency mission, will be used to characterize and analyze all possible 
sources of disturbance which perturb free-falling test masses from their geodesic motion. The system is composed of a 
single spacecraft (SC) enclosing a scientific payload, the LISA Technology Package (LTP), which is composed of two 
test masses (TMs) whose position is sensed by an interferometer. The spacecraft cannot simultaneously follow both 
masses, therefore the trajectory of only one test mass is used as a drag-free reference along the measurement axis. In 
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order to prevent the trajectories of the test-masses from diverging, the second test mass is capacitively actuated to 
follow the first (free-falling) TM. In the main science operating mode, the first interferometer channel measures the 
displacement of the SC relative to the free-falling TM. The second interferometer channel (the differential channel) 
measures the relative displacement between the two TMs. 

LPF is a controlled system which can only be fully assessed during flight operation, therefore a considerable number 
of experiments will be devoted to the identification of the details of the dynamics of the system. A dynamical model of 
LPF is built in advance on the basis of physical considerations and from the results of test campaigns. The dynamical 
model is parametric so that it can be updated on the basis of the experiments that will be conducted during mission 
operations. The overall aim of the process is to reach the best free-fall quality (below 3 x 10~^'' m s~^/\/Hz at 
frequencies around 1 mHz) in a step-by-step procedure in which the result of the previous experiment is used to 
adjust the detailed configuration of the following experiments [IHl]. 

Such a demanding program requires daily analysis of the instrument signals constrained by two major factors: 
i) the amount of available data is tightly constrained by LTP mission duration (90 days), the telemetry bandwidth, 
and the length of each data segment (typically 24 h); ii) the scientific interest is mainly focused on the analysis of 
noise sources which act directly on the TMs since that should provide a baseline reference for the forthcoming space- 
based gravitational waves observatories ^-[Sj . The direct forces on the TMs are expected to dominate the instrument 
output in the frequency range [0.1,10] mHz. Sample power spectra are typically calculated with Welch's averaging 
periodogram method (WOSA) f9^. In order to keep enough frequency resolution at low frequencies, the sample power 
spectra can not be averaged more than few times (we average 4 times in the present paper), this results in highly 
non-Gaussian data for which we are developing dedicated techniques. In particular the paper aims to propose a 
solution for two major data analysis challenges encountered in LPF: i) Different measurements of the same physical 
quantity can exhibit different noise content if they are performed under slightly different environmental conditions. 
The objective of LISA Pathfinder data analysis during operations will be to discover such differences, understand 
their origin and adjust spacecraft physical parameters accordingly. Such a problem requires reliable excess noise 
detection procedures which have to be based on solid statistical considerations; ii) Along with the demonstration of 
unprecedented test-mass free-fall, LPF will provide a model for the expected test mass force noise for future space- 
based gravitational wave detectors. In order to do this, wc need to be able to match an analytical model to a noisy 
power spectral density measurement. The quality of the match must be statistically quantified. Both data analysis 
problems deal with sample spectra and the corresponding statistical properties. Section |ll] reports on the properties 
of different experimental procedures for the detection of noise excess. In particular we considered two cases where 
the noise excess is evaluated with respect to reference data or a reference model. The accuracy of the methods is 
theoretically analyzed for the case of broadband and band-limited excess noise. In Section III the problem of noise 
parameter estimation for non-Gaussian data is explored and an algorithm based on maximum likelihood is derived. 
In parallel, the application of the Kolmogorov-Smirnov formalism to the construction of a goodness of fit estimator 
is discussed. Section |TV] reports briefly on the extension of the analysis procedures to the case of non-stationary noise 
and time- frequency investigations. In Section |V| we provide an application of the developed algorithms to synthetic 
data with LPF-like qualities. This allows us to shed light on the effects of data correlations on the accuracy of the 
developed excess noise estimators. The analysis procedures are available as MATLAB tools in the framework of the 
LTPDA Toolbox [TUHIS]; an object oriented MATLAB Toolbox for advanced data analysis. 



II. QUANTITATIVE DETECTION OF NOISE POWER VARIATIONS 

The problem of the detection of noise power variation in consecutive measurements can be formulated in two 
different ways, i) Two different measurements of the power are compared; ii) The different measurements of the 
power are compared with a reference model. The problem in the first case is of general character and can be applied 
to a wide range of experiments, the second case, on the other hand, assumes that a reference model for the noise power 
is known and data must be compared against the given model in order to understand if the system is performing 
under known conditions. The latter case is likely to be the scenario for LPF operations. The quantity typically used 
for the detection of noise power variations is the total energy content of the data series. It is defined as £ = Tj^i kil^- 
Unfortunately, £, provides a poor estimator for two reasons: i) As soon as the data series xq, . . . , xjy-i departs from 
zero mean Gaussian white noise, the statistic of £ becomes ill-defined and the definition of a confidence interval 
becomes cumbersome; ii) £ provides global information as it is not sensitive to noise changes in a given frequency 
band. While the first problem could be overcome (without little difficulty) by a numerical identification of the expected 
statistic, the second problem suggests that a spectral based estimator would provide supplementary information which 
could be fundamental in discriminating different noise sources. 
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A. Detection of noise variations with a model 



In the case that a reference model is available, the detection of excess noise in the spectral domain can be effectively 
implemented with a test on the normalized Welch's overlapped segment averaging (WOSA) spectrum i?woSA (/fc) 



defined in equation (A9). In particular, the test can pursue two different philosophies of which one aims to test a 
global scalar indicator of the properties of the data and another aims to test the details of the statistical distribution 
of the data. 

A sensible estimator for the first approach is provided by the integral of the normalized spectrum which, in the 
discrete case, can be written: 

Nf 



IR = ^i?WOSA(/fc). (1) 



fc=l 



In the simplifying assumption of independent spectral data, the statistic of each element of i?woSA i fk) is described 



by a gamma distribution as defined by equation ( A8 1 with S ~ I /Ns and h ^ Ng. The sum over the different values 
at the frequencies fk is still a gamma distribution with 5 = l/Ng and h = NgNf. The expectation value for IR is 
easily obtained as E [IR] = 5h = Nf. 

The natural estimator for the second approach is provided by the empirical cumulative distribution function 
(ECDF). The ECDF for a data series can be defined as F/v {x) = z/N, where z is the number of observations 
reporting a value less than or equal to x. x denotes the values taken by the data, in our case i^woSA jfk) with 



k = l,...,Nf. The ECDF can be tested against the theoretical expectation provided by equation (A8). If the 
model well represents the given sample spectrum then i?woSA (fk) is distributed according to the expected gamma 
distribution. Alternatively, if the sample spectrum contains a excess noise with respect to the reference model, then 



the distribution of the normalized WOSA spectrum will be different from the one reported in equation (A8). The 
hypothesis that the two distributions are equal can be tested if a 'distance' between the ECDF and the theoretical 
reference, Fr, is defined as: 

dK (x) = \Fm (x) - Fr {x)\ , (2) 

with cIk (x) assuming values on the interval [0, 1] and K ~ Nf. Kolmogorov found that the statistical properties of 

= max [dx (x)] (3) 

are independent from the specific distributions under test |13[ 114] . This property qualifies dx as an excellent 
candidate for the construction of a general test for cumulative distribution functions, the limiting statistic for dK was 
identified by Kolmogorov himself and then inverted by Smirnov [151 116j who provided an analytical expression for the 
calculation of dx as a function of the significance level. General details about the Kolmogorov-Smirnov test (KS-test) 
are reported in appendix [B| 

It is interesting to calculate the expected sensitivity for the two estimators. IR is expected to be distributed 
as a gamma distribution (at least when the model and the data are in agreement), the corresponding cumulative 
distribution function (CDF) is provided by the incomplete gamma function V {h, x/S) [T7]. The CDF assumes values 
in the interval [0, 1] as it defines the probability associated with the observation x. The inverse of the CDF provides 
the critical values, Xp, associated with a given probability, p e [0, 1]. The confidence range for a probability, p, for a 
gamma distributed variable can then be defined by the boundaries xi^ = 'P~^ (ck/2) and x^p = (1 — Q;/2), with 
a = 1 — p. We assert that the measured sample spectrum is compatible with the reference model if xi^ < IR < Xup 
for the given probability, p, or significance level, a. 

If the noise excess is provided by a scale factor, 7, which affects the noise on the complete band of frequencies, 
the expected value for IR changes to E (IR) — ^Nf . Therefore the detection threshold for 7 is fixed by the interval 
xiw/Nf < 7 < Xup/Nf. In other words, the IR can detect a noise difference with respect to the reference model 
only if 7 < xiw/Nf or if 7 > Xup/Nf. If 7 is non-zero only in a restricted band of frequencies [fa,fb], then the 
expectation value for IR changes to E (IR) = Nf — Nat + NatJ- In this case IR can detect the noise difference only 
if 7 < {xiw — Nf + Nab) /Nab or 7 > (xup — Nf + Nab) /Nab- Nab is the number of frequency points in the interval 
[fa, fb]- 

In the case of the Kolmogorov-Smirnov (KS) estimator, the ECDF of the WOSA normalized spectrum is compared 
with the theoretical expectation V {Ng,xNs). In the case that the noise excess is simply a constant scale factor, 7, over 
the whole frequency band, the expected CDF for the normalized WOSA spectrum is V {Ns,xNs/j). The expected 
value for the KS test variable is then written as: 



dx (7) = max 



V[Ng,^]~nNg,xNg) 



(4) 
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FIG. 1: Non-detection ranges for IR and KS estimators. Noise excess is assumed as a constant multiplicative factor 

7 extending along the whole band of frequencies. Darkened intervals define the thresholds for 7 detection, if 7 is 
larger or smaller than the shaded values it can be detected with a confidence of 95%. Ng refers to WOSA averages. 

Nf is the number of frequency points in the spectrum. 




FIG. 2: Non-detection ranges for IR and KS estimators in the case that noise excess coefficient 7 is different from 
zero in a restricted band [fa, fb]- Data are presented as a function of the ratio Nab/Nf where Nab is the number of 
frequency points in the interval [fa, fb] and Nf is the total number of frequency points in the spectrum. Ng refers to 
WOSA averages. Confidence level for excess detection is fixed to 95%. 

Once a significance value a is provided, the corresponding critical value, dK {ol), can be calculated from the equations 
for the inversion of the limiting CDF for (Ik [HI US] ■ The two distributions in equation Q are incompatible at the 
given significance level if (Ik (7) > dK (a). This defines the detection threshold for 7. If 7 is non-zero only in a 
restricted band of frequencies [fa,fb], then the expected distribution for the normalized WOSA estimator is difficult 
to calculate. Nevertheless the detection threshold for the KS estimator can be calculated numerically from synthetic 
data. 

In figure [Tl the ranges of non-detectability of the KS and IR estimators are reported as a function of Nf for three 
different values of the WOSA averages Ng. Data refer to the case that the scale factor, 7, extends over the complete 
frequency band. As can be clearly seen, the IR estimator always has a better sensitivity than the KS estimator. The 
sensitivity for the case of a band limited excess noise is reported in figure [2] as a function of the ratio Nab/Nf. As in 
the previous case, the IR estimator provides a better sensitivity with respect to the KS estimator. The difference is 
particularly relevant for Nab/Nf < 0.07, below such values the sensitivity of the KS estimator becomes poor. Both in 
figure [1] and in figure |2] the confidence level for 7 detection is fixed at 95%. 
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FIG. 3: Non detection ranges for the KS estimator for the comparison of two sample spectra. A) Noise excess 
coefficient 7 7^ on the full frequency band. Ng is the number of WOSA averages. Nf is the number of frequency 
points in the sample spectrum. B) 7 7^ only in a restricted interval of frequencies [fa, /&]■ Nat is the number of 

points in the interval. 



It is worth noting that KS algorithm can be used directly on time series to quantitatively gauging departures 
from a given distribution (e.g. Gaussian). Once the ECDF for the data is calculated, it can be compared with the 
expected distribution by calculating dxix) from equation [2] Since the distribution of dxix) coefficients is known it is 
straightforward to set a confidence threshold. The procedure can be applied even in presence of correlations thanks 
to the generalizations discussed in section |VB| and in appendix [B] 



B. Detection of noise variations without a model 



In the case that an excess of noise has to be detected by comparing different measurements the Parseval's theorem 
[5] (f = X^^o ^ (/fe)) suggests that the sum of the elements of a sample spectrum in a given frequency band 
Sab = Sfc=a ^ if'') could providc a sensitive estimator for noise power variations. Such an estimator would be loosely 
equivalent to the IR discussed above, except that its statistic is hard to determine in a typical experimental situation. 



The statistic of P (fk) at each frequency, fk, is defined by equation (A4). Therefore, in the case of non- white noise, 
its parameters depend on S (fk) and the statistic is different at different frequencies. Thus the statistic of Sab is not 
easily known. 

An interesting alternative to Sab is provided by the KS estimator. Given two data series {xn^} and {^n^} of length 
A'^i and N2 respectively. The hypothesis that their ECDFs have the same limiting cumulative distribution function 
F (x) can be tested if a distance in the ECDFs space is defined as: 

dK (x) ^ - FnAx)\ , (5) 

with dx (x) defined on the interval [0,1] and K — {N1N2) / {Ni + N2) [H]. Also in this case, the statistical 
properties of 

dx ~ max [dx (x)] (6) 

are independent from the distributions of Xn^ and . The same equations used in the case of the comparison with 
a given model can now be used for the inversion of the cumulative statistic of dx |3I]- Considering the simplifying 
assumption of independent spectral data, the sensitivity of the KS estimator at a given significance level can be 
calculated in analogy to what was discussed in the previous paragraph. 

In figure [3] the calculated interval for non-detection is reported for the case that the excess, 7, is extending over the 
whole frequency band and for the case that 7 7^ in [fa , fb] ■ 
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III. NOISE MODEL IDENTIFICATION 



Closely related to the problem of excess noise detection, the problem of noise model identification is one of the 
principal scientific objectives of LPF mission. Of particular interest is the identification of a model for the force noise 
acting on the TMs, which can be used as a baseline for the forthcoming space-based gravitational wave observatories. 

The main constrains on the identification of force noise on the TMs in LPF are fixed by i) the limited data series 
length, typically 24 h; ii) the frequency range in which force noise on the TMs is dominating the signal is below 
10 mHz; in) force noise data are not directly accessible since the system measures and reports TM displacement, 
requiring that force noise on TMs be reconstructed by a numerical procedure [THl US] • 

The result of the combination of the first two constraints is that the number of segment averages in the WOSA 
procedure for sample spectrum estimation should be taken as low as possible so as to have a reasonable number 
of frequency points in the range / € [0.1, 10] mHz. As a consequence, the distribution of the WOSA spectrum 
strongly departs from a Gaussian distribution |32) . meaning that the classical least-squares minimization procedure 
for parameter estimation is not well conditioned and a full maximum likelihood procedure is required. 



A. Maximum Likelihood Parameter Estimation 



If we replace S (fk) in equation ( A9 ) with a parametric model for the spectrum, S [fk] 9i, . . . , Oh), the normalized 
WOSA spectrum becomes parametric, RwoSA {fk]di, ... ,6h), and can be used for the estimation of noise model 
parameters {^i, . . . , Oh}. In this section we develop the likelihood formalism for the simple case that the noise model 
is a function of a single parameter, S {fh]9) — 0Sq (fk)- This allows us, not only to find a sensible goodness of fit 
estimator for i?wosA {fk',0) that can be used also in the case of multiple parameters, but also to place the excess 
noise estimator, IR, in a more solid theoretical framework. 

Indicating with ^true the 'true' value for the parameter, i?woSA (A-; &) can be rewritten as: 

n / . ,x , AvoSA ifk) 

-KwosA [jk;(p) = <p 



^true-So (fk) (7) 

= 0i?WOSA ifk) ■ 

Here (j) = ^truo/^- Assuming that 6'true 5*0 (/fc) correctly reproduces the expected value for the spectrum, the 
distribution of i?\voSA ifk) is reported in equation (A8) with 5 = l/Ng and h = Ng. The distribution of the samples 
■RwoSA ifk; 0) is then: 

Under the simplifying assumption that the values of -RwoSA ifk', <t>) are independent for different fk, the likelihood 
function can be written as: 

C{h,5,ci,)^\{f{yk;h,5,4>)l^y. (9) 

k 

Here Ay is a constant term required to have a finite probability from the probability distribution function / (y^; h, 5, cp). 
Uk are observed samples corresponding to i?wosA (/fc! '/')■ It is typically more convenient to work with the natural 
logarithm of the likelihood function I {h, S, (p) = ln£ {h, S, (p). 

I {h, S, (t>)^{h-l)J2 In Vk-^Y. -NfhXn cp. (10) 

fe k 

Nf is the total number of frequency samples. Taking the first derivative with respect to (p and equating to we find 
the maximum likelihood estimator for the parameter (p: 

A= ^^i?wosA(A-;0)-iV/. (11) 

k 

The value of (p at which A = corresponds to the maximum likelihood estimation for the parameter. The A 
estimator is unbiased, since, remembering that (p = Otmc/O, it can be verified that Vmig^e^^^^ E[K\ = 0. For the 
practical purpose we can find the zero crossing of the reduced estimator 

A = ^i?wosA(A;e)-iV/, (12) 
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which is crossing zero at the same 9 value of A, since (/>—?> 1 when 9 — >■ 0true- It is worth noting that A is practically 
the IR excess noise estimator with the expectation value Nf subtracted. It is then worth noting that A can be used 
not only in the simple case of one parameter but it can also be applied in the general case of a model S {fk;9i, . . . ,9h) 
since the condition A — >■ when {^i, . . . , 9h} — > {^itruc, • • ■ , ^//truc} is always satisfied. 



B. Kolmogorov-Smirnov Parameter estimation 



As discussed in section |II A[ the KS estimator can be used as an alternative to a maximum likelihood procedure 
for parameter estimation. 

Thanks to the statistical properties of i?wosA jfk', 9), the closer 9 is to 6'truc, the better the distribution of 



-RwoSA (/fc! 6') is described by equation (A8) with d = l/Ng and h — Ng. Therefore, the Kolmogorov-Smirnov 



distance parameter provides an effective goodness-of-fit estimator: 

dKi9)^ma^\FRix;9)-riN,,xN,)\. (13) 

Fa (a) is the ECDF for the current i?wosA ifk', 9) estimate, V {Ng, xNs) is the limiting distribution for 9 — > 6'truo- 

KS estimation for 9 is obtained by the minimization of (9) with respect to 9. A confidence range for the 
parameter estimation can be readily defined from the non-rejection region of the KS-test at a given significance level. 
In practice, having defined a significance level, the corresponding critical values of the KS statistic, dx (a), can be 
calculated with standard equations [TB] or Monte Carlo simulations in the case of correlated data. The values of 9 
for which dK {9) = dx (a) provide the boundary for the confidence range at the given significance. The KS statistic 
can also be successfully applied to multi-parameter identification, since the convergence of Fj^ {a) to V {Ns,xNs) is 
always verified when {6*1, ... , 9h} {6'itruo, • • ■ , ^'fftruc}- 



IV. ANALYSIS OF NON-STATIONARY NOISE 



The implementation of noise analysis procedures was so far discussed in the context of stationary or pseudo- 
stationary noise [33] . In the case of truly non-stationary noise the spectral content of a time series is investigated 
by time-frequency analysis techniques which include the spectrogram and the wavelet transform. The spectrogram is 
estimated by the square modulus of the short-time Fourier transform of the data |20) . It provides a direct extension 
of the PSD formalism to non-stationary time series. Given a data series of N samples, a fraction of length Q < is 
windowed and the Fourier transform is applied. Then the data window is time shifted and the process is repeated. The 
calculation of the spectrogram is based on data windowing and the application of the Fourier transform, therefore the 
considerations noted in the previous sections for stationary noise can be applied directly to the spectrogram analysis 
of non-stationary noise. 

Since the short-time Fourier transform has the same resolution across the time- frequency plane, it is often preferable 
to resort to the wavelet transform. Wavelet transform is a decomposition of the time series over time-frequency 
elements that are obtained by scaling and translating a mother function G (M): 



^u,s = . (14) 

The function tp (t) has zero average and the wavelet elements are normalized to 1. The wavelet transform of a 
function f{t) is then defined as: 



Wf{u, s) = p fit) (^-^)dt. (15) 

In the discrete case, the result of the wavelet transform on a time series is an array of coefficients Wu,s where u is 
the time index and s is the scale index which is associated to a given frequency band [H] . In the case of uncorrelated 
Gaussian noise, the distribution of the coefficients Wu,s is still Gaussian with a certain amount of correlation introduced 
by the convolution-like transform operation |21j . In such a favorable situation the extension of the method described in 
the stationary case it appears straightforward. In particular Kolmogorov-Smirnov procedures are excellent candidates 
since their robustness to correlation and the possibility of extending KS distance definition to a two-dimensional space 
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Frequency [Hz] 

FIG. 4: Projection of the spectrum of Qcs- 'Readout' is the projection of the readout noise to Ooff . 'TMs' is the 
projection of the force noise on the test masses to Ooff . 'TMs + Readout' is the complete noise projection for Ooff , it 
represents the basehne noise level assumed in the present paper. 'LPF worst case' refers to a worst case scenario for 

Oeff and 'LPF Spec' corresponds to the mission specifications. 



V. APPLICATION TO LPF SYNTHETIC DATA 

In this section the different procedures for excess noise detection and parameter estimation are applied to synthetic 
LPF data. This provides not only an interesting framework for testing their accuracy and precision, but it also helps 
clarifying the role of correlations among spectral data. 

A. Synthetic data and noise projection 

LPF provides two output channels along the principal measurement axis which are sensing the displacement of the 
SC relative to the free falling TM and the relative displacement between the TMs. 

From the knowledge of the displacement signals and a linear model for the system dynamics, an effective force-per- 
unit-mass, Oeff, acting on the TMs can be extracted. Oeff is the combination of the 'true' force-per-unit-mass acting 
on TMs and a projected interferometer readout noise. Details of the calculation are reported in appendix [P} 

Following the same scheme, a given prediction for the different input noise sources can be projected to a prediction 
for the power spectrum of Ocfr, as reported in figure [4] Here the model for the spectrum of force noise acting on TMs 
and the model for the spectrum of readout noise are projected into a model for the power spectrum of Ocff (TMs + 
Readout in the figure notation), which represents this paper's baseline for the spectrum of Ucb- 

In figure |4] we also report the project specifications for LPF and the expected noise spectrum for Ooff in a worst 
case scenario. In our baseline we assumed a reduced force noise on the TMs compared to the worst case but choose 
to keep the worst case for the readout noise. This was done in order to represent one of the possible scenarios (not 
the best one) that can be experienced during the mission. 

The model assumed for the force noise on the TMs is characterized by a low frequency 1//^ behavior and a flat part 
for f > 1 mHz. The model can be written as ^tm (/) = (^S^tm (/)■ ^ is an adjustable parameter which assumes values 
= 1 for the worst case scenario and 6 = 0.1 for our baseline model. It is worth noting that Stm (/) is projected 
(together with the readout noise model) through LPF dynamics in order to obtain the expected noise spectrum for 
Ooff which we indicate with Sa (/). Sa (/) with 6 = 0.1 corresponds to 'TMs + Readout' in figure|4j Sa (/) with 6 = 1 
corresponds, instead, to the 'LPF worst case'. 

B. Excess noise detection 

A change of the noise level on S'tm (/) {(^ 7^ 0.1) produces a variation of the energy content of Ocff . Such variation, 
which may be 'improperly' identified as excess noise, can be detected with the procedures defined in section [iTj In 
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particular we tested the detection of excess noise between two data series and between a data series and a reference 
model. Synthetic data were produced according to the following procedure: 

1. Different models for Sa (/) are produced changing 9 around the reference value 6 — 0.1. Readout noise level is 
kept fixed. 

2. Corresponding noise time series for a^g are generated using the procedure reported in |24j . The time series are 
24 h long and have a sampling frequency of 1 Hz. 

3. Sample spectra are calculated for each series with the WOSA algorithm. We chose a Blackman-Harris data 
window, 50% segment overlap and number of segments averages Ng = 4. 

4. The analysis is restricted to the frequency interval [0.1, 10] mHz since, as can be seen from figure |4j it represents 
the region in which the force noise on the TMs dominates Sa (/). 

Spectral data are tested for excess noise. We used the KS algorithm (equation (|6])) in the case of the test of two 
data series. The data series for 9 ^ 0.1 are compared against the reference series with 9 = 0.1. In the case of the test 
of a data series against a model, both the KS algorithm (equation ([3|) and the IR algorithm (equation Q) are used. 
The reference model is the projected Sa (/) for 6* = 0.1. The results of the tests are summarized in table urKS critical 
values dx (a) and IR confidence intervals are calculated for a significance level a = 0.05 which corresponds to a 95% 
confidence. 

Each value of 9 corresponds to a value of the in-band energy content E (9) of Sa (/) in the analyzed frequency band. 
We report in table|l]the relative change in energy AE/E corresponding to a relative change in 9. The quantity AE/E 
plays the same role of the parameter 7 in figures [T] and [3j even though the two quantities are not completely comparable 
since 7 in figure [T] assumes independence of the data. Spectral data are correlated among different frequency values 
because of two effects [9l [25] : i) Data windowing which corresponds to a convolution in the frequency domain of the 
window function with the sample spectrum; ii) WOSA overlapped segment averaging. The first effect is unavoidable, 
the second effect, instead, can be attenuated by a proper choice of the segment overlap. It can be demonstrated [M] 
that for a Blackmann-Harris window the effect is negligible with 50% overlap. 



TABLE I: Detection of noise differences in the frequency band [0.1, 10] mHz. The symbol y indicates compatibility 
between tested objects. The symbol x is instead used for indicating test rejection. 9 is an adjustable parameter 
which assumes value 9 = 0.1 for our baseline model. Different values of 9 correspond to different values of the 
in-band energy content E (9) of Sa (/). We reported here relative values with respect to reference. Details on 
^'k (Q^)' "^^^ i^) MC confidence interval for IR estimator can be found in appendix fs^ and Icj respectively. 
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Since the standard statistics for the estimators {dx (a) for KS and confidence interval for IR) are calculated in 
section |ll with the assumption of data independence, the standard critical values, dx (a), for the KS estimators and 
the confidence intervals for the IR estimator can be applied only if the correlations among data are negligible. If this 
is not the case, the effective statistic of the estimators can be numerically calculated with JVIonte Carlo simulations. 
The corresponding results of a Monte Carlo simulation with A^mc = 5000 realizations of the reference data series are 
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indicated in table|l]with the sufRx MC. The symbol ^ is used to indicate that the spectral data for the corresponding 
i^E / E is compatible with the reference. On the contrary the symbol x indicates a rejection. 

Observing the test results reported in columns 7 and 9 for the KS algorithm and the results in columns 11 and 
12 for the IR algorithm it is readily seen that correlations among data play a role. In the case of the KS test, the 
comparison of the data with dK (a) determines a rejection for /S.E/E — 0.05. On the other hand, the same value is 
accepted when Monte Carlo result {a) is used. In the case of the IR estimator, the comparison with the standard 
confidence interval leads to the rejection in correspondence with /S.E/E = 0.05 and l^E/E = —0.04. Such values are 
instead considered compatible by the Monte Carlo confidence interval. These results provide us with clear information 
that the presence of correlations among data has affected the tests statistics and therefore the standard equations, 
assuming independence among data, are not usable in this situation. 

Looking at the results for the KS test between two data series, we discover that dK (a) and (a) are practically 
equal. In fact the results for the two corresponding columns of table |l] (columns 4 and 5) are in perfect agreement. 
This is the practical result of one of the most interesting properties of the KS test. The test is based on equation (|6|. 
It states that the two empirical cumulative distributions under test have the same limiting CDF. Since the sample 
spectra in our test are calculated following the same WOSA procedure, the degree of correlation among different 
frequency points is the same, therefore the test statistic is not spoiled. 

Comparing dK (a) and d^'~^ (a) of columns 7 and 9 we see that the effect of correlation is to increase the maximum 
expected spread between ECDF and limiting CDF. Therefore the effect of data correlation is to distort the expected 
statistic for dK- The 'distortion' of dK statistic can be taken into account if an effective value for the parameter K is 
introduced. In the case of the comparison of the ECDF for correlated data against a theoretical CDF the application of 
the standard values for dK, where K = Nf, Nf being the number of data elements, leads to a statistically unfair test. 
We then discovered that test fairness can be recovered if an effective value for K is used rather than the standard 
K = Nf. In particular, for spectral data produced with the WOSA method, Blackmann-Harris window, A^^ = 4 
averages on 50% overlapped segments, we obtained KeS = (3Nf with P = 0.55 for a significance level a = 0.05. It is 
worth noting that the value of (3 is independent from the number of data points considered, and from the spectral 
shape, provided that the different shapes have reasonably comparable smoothness on a frequency interval comparable 
with the width of the first lobe of the data window. As an example, the value of /3 = 0.55 is valid for LPF-like data 
and for white-noise equivalently. The requirement on the smoothness of the spectrum is connected to the expression 
of the correlations introduced by data windowing. It can be demonstrated that if the spectrum can be assumed 
constant in a frequency range of the order of the width of the first lobe of data window [^9 then the correlations are 
independent from the particular shape of the spectrum and are determined only by the window function. For such 
class of spectra we expect the same value for /? once the required significance level is fixed. 

C. Noise model identification 

The problem of noise parameter identification is strictly connected to the problem of excess noise detection by a 
comparison of a data series with a reference model. While in excess noise detection different data series are compared 
with a given reference, in parameter estimation different realizations of a parametric model are compared with a 
dataset in order to find the best fit. Due to this, the same algorithms (i.e. KS and IR) can be applied to the solution 
of the two problems. Precision and accuracy in the estimation of the parameter, 9, controlling the excess force noise 
on the TMs is tested with a Monte Carlo simulation over A^mc = 5 x 10'^ realizations of the same process. The data 
series reproduce Ucs corresponding to the reference value 9 — 0.1. Sample spectra are calculated with the procedure 
described above and the analysis is restricted to the frequency range [0.1, 10] mHz. For each realization, data are 
compared with a bank of models obtained by the projection of TMs force noise ^tm (/; 6) = ^"^tm (/) ^'^d readout 
noise into Sa (/; 9) for different values of 9 around the reference value. The KS and IR estimators are calculated 
for each element of the model bank, in particular it is simpler to analyze the results in terms of IR = |IR— A/|. 
Both KS and IR are expected to have a minimum corresponding to the best estimate for the parameter 9. The 
two methods proposed here are compared with the performance of a classical weighted least-squares method, which 
works by minimizing the mean squared error MSE — J^f,, {{PwoSA (fk) — Sa {fk',9)) /Sa {fk',S))'^, and a Hubcr's 
norm estimator (details are reported in appendix |e| . The results of the analysis are reported in figure [s] where the 
histograms of the best fit vales over A'mc realizations are reported for the four procedures. In the same figure we also 
show the evolution along the model grid of the four estimators for a particular set of data from the available A'mc ■ 

The distributions for the best fit parameter are reasonably symmetric for all the estimators; mean values and 
sample standard deviations are respectively 9ks = 0.100, aKS — 0.012, Oir — 0.100, am = 0.012, 9mse = 0.148, 
o'MSE = 0.016, Onuber — 0.130 and cjHuber = 0.014. From the analysis of the Monte Carlo results, it is readily seen 
that both the KS and the IR estimators provide equivalently precise and accurate results. On the other hand, the 
MSE algorithm provides a poor estimation, both from the accuracy and from the precision point of view. The best 
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FIG. 5: KS, IR, MSB and Huber's norm accuracy and precision in the estimation of the force-pcr-unit-of-mass noise 
parameter 0. Corresponding histograms report the result of -/Vmc = 5000 MC reahzations of the parameter search on 
a grid of template models. We also show (black dots) the evolution along the model grid of the four estimators for a 

particular set of data taken from the A'mc available. 



estimation for the parameter is Omse = 0.148, which is strongly biased with respect to the reference value oi 9 = 0.1. 
Also, the distribution of the parameter values is wider than those obtained from the KS and IR estimators. Huber's 
norm estimator, with the chosen parameter c = 0.05 |Ej performs better than MSE but the result 0Huber — 0.130 is 
still far from the true value. 

It is worth discussing the effect of correlations among spectral data on the estimation of the parameter Q. As can 
be observed from the results of the MC simulation, the accuracy of the estimators is not affected by correlations; 
the results for KS and IR estimators are practically indistinguishable. Some problems with IR can arise from data 
correlations when the confidence interval for a single estimation is required. As discussed in the previous paragraph, 
correlations modify the statistic of the IR estimator making it impossible to easily calculate confidence intervals from 
the standard equations. The statistic of the KS estimator is affected by correlations too, but we have demonstrated 
that accurate critical values (a) can be recovered if an effective value for K is used, ifos — PNf. /3 is a shape 
parameter depending only on the correlations in the spectral data and on the required significance level. It can be 
calculated for a specific spectrum (e.g., white noise) and effectively extended to a wide family of spectral shapes. 
Using the value /3 = 0.65 reported in appendix |Bj we obtain dP^ {a) ~ 0.064 for a = 1 — 0.68%. The intersection of 
such a value with the curve of dx as a function of 9 reported as black dots in figure [5|\ provides a 68% confidence 
interval for the single estimation. In this specific case such an interval is 6* € [0.086,0.114]. 



VI. CONCLUSIONS 



The problem of excess noise detection and noise parameter estimation for non-Gaussian data is analyzed in the 
framework of the LISA Pathfinder mission. Excess noise detection can be approached in two ways. In one way, 
the noise content of a data series is compared with a reference data series, in the other way the noise content of 
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the data is compared with a reference model. In the first case, simple estimators like the total energy content in a 
data scries are not suitable for formulating quantitative statements on a solid statistical basis. As an alternative a 
Kolmogorov-Smirnov (KS) estimator is proposed and successfully applied to LPF synthetic data. The KS estimator 
has the advantage of being independent of the statistical properties of the data under test. It is demonstrated in 
the paper that such a convenient property makes the estimator robust to correlations among spectral data. Two 
different estimators arc investigated for the problem of comparing a data series to a model for excess noise detection. 
One estimator (IR) is based on the integral of the normalized WOSA spectrum, the other is a KS estimator for the 
comparison of an empirical cumulative distribution function with a limiting theoretical function. Despite the fact 
that IR estimator proves to have better sensitivity on independent data, the versatility of the KS estimator is highly 
advantageous on correlated data. Since the statistics of the estimators are distorted by the presence of correlations 
among test data, the standard and simple procedures for the determination of the confidence intervals, based on the 
inversion of the limiting distribution function, provide inaccurate results. While, in the case of the IR estimator, the 
problem can be overcome only with dedicated Monte Carlo simulations, we demonstrated that the introduction of 
a shape parameter allows us to use standard equations to calculate proper boundaries for KS confidence intervals. 
The shape parameter depends on the required significance level and on the data correlations. It does not depend 
on the number of data considered for the test. Since correlations among spectral data are mainly introduced by the 
windowing process, the shape parameter is fixed for a wide family of spectral shapes. For example, synthetic LPF 
data share the same shape parameter with white noise. 

Closely related to noise excess detection, the problem of noise parameter identification is analyzed with a maximum 
likelihood approach which, in the particular case of linear dependence on a single parameter, provides an algorithm 
analogous to the IR estimator used for excess noise detection. A KS algorithm was proposed as an alternative to the 
IR algorithm and the accuracy and precision of both were tested with a Monte Carlo simulation on LPF synthetic 
data. Both the IR and KS estimators were demonstrated to give equivalently good results, even though the capability 
of the KS to handle data correlation is a clear advantage for the definition of a confidence interval for the estimated 
noise parameter. 

Data analysis procedures introduced in this paper are easily extended to the vast context of time-frequency analysis 
of non-stationary noise. KS algorithm can be applied effectively both to spectrogram and wavelets coefficients provided 
that correlations among data are taken into account. KS algorithm has the advantage that can be generalized to two 
dimensions thus allowing to extend the analysis to the time-frequency plane. 

Appendix A: Statistic of the sample spectrum 
In the case of a discrete, real-valued stationary process {xh}, the continuous spectral density function is defined as: 

oo 

S{f)=T ^ s^e-'2.AT_ ^Al) 

h=—oo 

Here Sh is the autocovariance sequence of the process {xh}, T is the sampling period and / is the frequency expressed 
in Hz. / is defined in the range |/| < Jnq = 1/2T and Jnq is known as the Nyquist frequency. In the case of a 
finite representation xq, . . . ,xn-i of the discrete process {xh.}, the approximation to the spectral density function is 
provided by the sample spectrum 

2 

(A2) 

/ in this case is also defined on the interval [— /wq, /wg]- If the sample spectrum is calculated on the grid of Fourier 
frequencies (/^ = k/ (NT), \k\ < N/2) then it corresponds to the squared modulus of the discrete Fourier transform 
of the data sequence xq, ■ ■ ■ ,xn-i- In practical applications only the positive frequency part of the spectrum is 
considered, and the one sided sample spectrum is defined as P (fk) = 2P (fk) with k = 0,1, . . . , N/2. In the rest of 
the paper the one sided sample spectrum will be simply named the sample spectrum. 

If the data series xq, . . . ,Xn-i is Gaussian distributed and the elements xj are independent then the Fourier trans- 
form produces a complex series X (/j.) whose elements are approximately independent and their real and imaginary 
parts are Gaussian distributed. The term \X (/fc)P = |3? [X + |Ss [X {fk)]\^ is then the sum of two independent 

variables distributed as a xt where u is the number of degrees-of-freedom of the distribution = 1 in this case). If 
the correlations among the elements of the data series xo,. . . ,xn-i are non- vanishing, the statistical properties of 
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the sample spectrum, in the simphfying assumption of independent P (fk) elements, can be calculated from the case 
of a xt distributed variable y multiplied by a constant z = Xy. The characteristic function for z is 

<j),{t) = E[e'*^y] = {l-2it\)-^ . (A3) 

Here E [] indicates the expected value. The inverse Fourier transform of 0^ (t) provides the probability density function 
for z: 

^ {^At)]^-—^--V-v- (A4) 



(2A)^ r I 



2y 



This is a gamma distribution, / (z; k, 9), with k = t//2 and 6 = 2A. F (fc) is the gamma function. In the case of the 
sample spectrum at a given frequency, z — P (fk) and X = E[P (fk)] /v ~ S (/) /v [35]. It is useful, for the statistical 
analysis of the spectrum, to introduce the normalized sample spectrum 

«(A)^.ff, (A5) 

which, at each frequency fk, is distributed as x^. 



The sample spectrum in equation (A2) can also be seen as a special case of the power spectral density, S (/), when 
the infinite data series {xh\ is chopped by a square data window of length N . This operation introduces a considerable 
amount of spectral leakage because of the convolution with the frequency response of the square data window [HI [2S] ■ 
Therefore it is common practice to multiply the time series a;o, . . . , xm-i with a more pcrformant data window, which 
increases the accuracy of the sample spectrum in the case of processes with a high dynamic range. The application of 
a data window introduces correlations among different elements of the sample spectrum. Such correlations affect the 
statistics of the spectrum, resulting in a change in the probability distribution of the sample spectrum. An analytical 
treatment of the spectrum statistics under such conditions is cumbersome, and it is easier to numerically evaluate 
(with a Monte Carlo simulation) the statistics of the sample spectrum for the case of interest. 

In order to improve the variance properties of the sample spectrum, Welch's overlapped segment averaging (WOSA) 
method is applied Jj. The data series a;o, . . . , xjq-i is divided in overlapping windowed segments. The estimates of the 
sample spectrum of each segment are then averaged. The practice of averaging overlapping segments can modify the 
expected statistics of the spectrum since the data in different segments can be correlated. Then, even in the simplifying 
assumption of vanishing spectral correlations from windowing and overlapping, the averaging process changes the 
statistics of the estimated sample spectrum. In the assumption of vanishing window and overlap correlations, the 
statistic of the WOSA spectrum corresponds to the average of Ng gamma distributed variables 

1 

PwosA(/fc) = ^E^j(^)- (A6) 

" J=l 

The critical function for the sum is 0p„osa (^) = Yij i^)' where (t) is the critical function for Pj/Ng. Thanks 



to equation ( A3 1 



'Apwosa it) = (^1 - ]^ j ' (A7) 

where 9 — 2S (fk) jv and k — 1^/2. The inverse Fourier transform of (j)PvjosA (^) provides the probability density 
function for the WOSA spectrum 

/wosA {z; h, S) = , (A8) 

with S = S (fk) /Ns and h = Ng. Again it is useful to define a normalized WOSA spectrum as 

p (f ^_ fwosA ifk) . . „^ 

^ {Jk) 



which is gamma distributed (equation (A8|) with S — l/Ng and h — Ng- 
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Appendix B: Kolmogorov - Smirnov Test 

Kolmogorov - Smirnov is a well known test for distributions [TSHIBl [27l [28] . An empirical cumulative distribution 
(ECDF) is tested against a continuous theoretical model or, alternatively, two ECDFs are tested with the hypothesis 
that they share the same limiting cumulative distribution function. Indicating with / (x) the probability density 
function associated with a given random process X, the corresponding cumulative distribution function (CDF) is 
defined as: 

F (x) = Prob [X <x]^ f f (u) du. (Bl) 

J —oo 

Given a particular realization of the random process X: 

Xn ^ {xi,. . . ,xn} , (B2) 

ECDF is written as Fj^ (x) — k/N where k is the number of observations which is smaller or equal to x. 

Given two data series Xjsi and Ym, with N and hi not necessarily equal, we can test if the two series are two 
particular realizations of the same random variable by analysis of their ECDFs. Under the hypothesis that the two data 
series comes from the same distribution function, Kolmogorov has demonstrated that the maximum distance between 
the two ECDFs has a limiting distribution which is independent from the statistical properties of the corresponding 
random variable. If the test is performed against a theoretical distribution, the distance is defined as: 

dK = max \Fn {x) ~ F {x)\ . (B3) 

In such a case K = N . In alternative, If the test is performed between two ECDFs, K — (NM) / {N + M) and: 



dK = max \Fm (x) - Fm [x) \ . (B4) 

The test is defined as follows: 

1. In the case of the test on a single data series, the null hypothesis is that the data are realizations of a random 
variable which is distributed according to the given probability distribution. In the case of two data series, 
the null hypothesis is that the two data series are realizations of the same random variable, which means their 
ECDFs should tend to the same CDF. The test rejects or accepts the null hypothesis on the basis of the analysis 

of djijyi. 

2. A significance level a is defined as the probability that the test rejects the null hypothesis when it is indeed 
true. 

3. The test can be formulated in terms of critical values. The critical value for the test is the value of dK (a) 
corresponding to the significance level. Then if dK > dK (a), the null hypothesis is rejected. 

KS theory was formulated for independent data sets and the available equations for critical values are valid only if 
this condition is satisfied [TB] . The test can be formulated also in the presence of data correlation but the distortion 
to dK statistic introduced should be taken into account. This is possible if an effective value for K is introduced 
as Kcs = l3K, with /3 a shape parameter depending only on the data correlations and the required significance 
level. Alternatively, realistic critical values can be calculated with dedicated Monte Carlo (MC) simulations [29]. 
The advantage of the shape parameter j3 is that it depends only on correlations and the required significance level. 
Therefore it can be determined for a specific spectrum (e.g., white noise) and shared among a wide family of spectral 
shapes. Once /3 is known it can be used to calculate critical values for correlated data using standard equations 
reported in the literature [16]. 

Focusing on the particular problem, we performed a Monte Carlo estimation of dK (a) for WOSA spectra repre- 
sentative of LPF. The number of frequency data considered is Nf — 341, corresponding K values are K = Nf in the 
case of the test against a theoretical distribution and K — Nf/2 in the case of the test between two ECDFs. The 
results are summarized in table jll] In the same table we report the values of /3 that are required to obtain proper 
critical values from the standard equations in the case of the test against a theoretical distribution. 
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TABLE II: Table of KS critical values for correlated spectral data. Critical values are calculated with a Monte Carlo 

simulation for different values of the significance level a. Datasets used are representative of the data analyzed in 
the current paper. We reported the values for testing an ECDF against a theoretical CDF [K = Nf) and the values 
for testing two ECDFs for the same limiting CDF (K = Nf/2). We also reported the values of the shape parameter 
beta that can be used to calculate proper critical values from standard equations in the case of correlated data, a 
refers to the significance level whereas 1 — a is the corresponding confidence level for the test. Nf = 341. 
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Appendix C: IR Test 

The statistics of the IR excess noise estimator are numerically estimated with a Monte Carlo simulation on the 
spectral data used for the present paper. Data series are 24 h long synthetic reconstructions of the force-per-unit-mass 
expected on the TMs. WOSA spectra are calculated with Ng — 4 averages on 50% overlapped segments which were 
windowed with the Blackman-Harris window. Results are reported in table |III| 

TABLE III: Tabic of confidence bounds for the IR estimator on correlated spectral data. The values are calculated 
with a Monte Carlo simulation for different values of the significance level a. Data sets used are representative of 
the data analyzed in the current paper. The number of frequency points is Nf — 341, corresponding to the number 
of available spectral data in the range [0.1, 10] mHz, for a time series 24 h long and number of averages for the 

WOSA estimator N^ = 4. 
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0.32 


0.68 


325.88 


357.41 


0.10 


0.90 


316.33 


367.28 


0.05 


0.95 


311.33 


372.55 


0.01 


0.99 


300.48 


380.69 



Appendix D: Conversion of displacement noise 

LPF can be considered as a three body controlled system composed of the two TMs and the spacecraft (SC). Its 
equations of motion along the measurement axis can be written as: 

mix'i + mix'sc + miLolxi — fi 

1712X2 + m2X"sc + m2i^lx2 = + fc2 (Dl) 
mscX'sc - miUjlxi - m2Ujjx2 = fsc - fc2 + fcsc ■ 

Here: 

• xi and X2 are TMs coordinates along the sensitive axis. They are relative coordinates in the SC reference frame. 

• Xsc is the absolute SC coordinate along the sensitive axis. 

• nil, fTi2 and vrisc are the masses of the two TMs and of the SC. 

• and UJ2 are the parasitic stiffnesses coupling the TMs and the SC. The TMs are coupled to the spacecraft 
through the parasitic stiffness thus producing an oscillator like equation of motion. The spacecraft at the same 
time experiences reaction forces given by —miuifxi and —m2Uj2X2- 

• fiT I2 and fsc are the forces acting on TMs and SC respectively. 
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• fc2 and fcsc are control forces on the second TM and the SC respectively. Since fc2 is an internal force to the 
system the SC experiences a reaction force —fc2- 

• Dots over symbols represent time derivatives. 

In the main LPF science operation mode, one TM (indicated here as TMi) is in free-fall and provides the reference 
for the other TM (TM2) and the SC. In order to avoid unwanted drifting, both TM2 and the SC are controlled to 
follow TMi. It is worth noting that the system is, by construction, symmetric and the role of the two TMs can be 
inverted. In order to avoid confusion, we indicate with TMi the free-fall reference and with TM2 the actuated TM. 

Moving to the Laplace domain, substituting for the SC dynamics and substituting for the differential coordinate 



XA of TM2 with respect to TMi, the equation (Dl) can be rewritten as: 



2 2/1 "^1 \ 2 ™2 2 ™2 

s XiOJ-, 1 H ] xi + ^1 + ^2 = 

\ rrisc J nisc rUsc 

fl fsC I 1 TJ ^ TT 

= + H2OA HscOl 

mi lUsc rUsc rrisc (Ij2) 

S X\ 4- {u)2 — UJi) Xi + L02Xl^ — 

f2 /l I 1 rr 

= \ -H20A- 

TO2 mi 7712 



Here: 



• oi and OA are output displacement signals as provided by the interferometer readout system, oi is the displace- 
ment between the SC and the TMi. oa is the displacement of TM2 relative to TMi. 

• H2 and Hsc are transfer functions of the control systems on TM2 and SC. The force applied by the controllers 
is calculated on the basis of the output displacement, therefore fc2 — H2OA and fcsc = HgcOi- 

Calculations can be more easily performed if we introduce a matrix notation: 

(D3) 



(D4) 






(D5) 



.... . ..J 



G = ( ° I (D7) 



C = ™- ™f c (D8) 



(D9) 
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With the notation introduced, the equation of motion can be rewritten: 

Dx = Gf + CHo. 



(DIO) 



Then we should consider that the output displacement o corresponds to the measurement of x provided by the 
sensing system: 



(Dll) 



Here S is a 2 x 2 sensing matrix and Om is the readout noise. 



Substituting equation (Dll) in equation (DIO I, the dynamics can be rewritten in terms of the known output o; 

)r„. (D12) 



D S 



o C H 



G f + D S 



The quantity G • f represents the force-per-unit-mass acting on the test masses. Such a quantity is not directly 
known from the system output since the available signal is o. From equation (D12) it is readily seen that an effective 



force-per-unit-mass acting on the TMs can be reconstructed from the knowledge of o if the force applied by the control 
system (C • H • o) is calculated and subtracted from the reconstructed dynamics D • • o. Therefore the available 
quantity is g = G • f -I- D • S^^ • Om which contains the force-per-unit-mass acting on the TMs corrupted by the 
readout noise of the system. 



The expected values for g can be obtained thanks to the dynamics reported in equation (D12) once the expected 
values for f and Om are known. The same applies to the projection of the noise spectra for f and Om to g. 



Appendix E: Huber's norm 

Ruber's norm [3U] is a way to construct a goodness of fit estimator which is more robust to outliers and non- 
Gausianity than the standard mean squared error (MSB). The norm is constructed as '^iPi^i), where are the 
residuals between a data series and a parametric model. The function p{ri) s defined as: 



p{n) = \^'"^ fo^l^*l<c (El) 

\c\ri \ - \(? for \ri\ > c. 

The value of the threshold constant c may depend on the given dataset, therefore some efforts must be 
spent to select the value of c providing the most accurate results. Normalized residuals are defined as Vi = 



{PwosA ifi) ~ Sa iff, 0)) /Sa [fi] 0) in accordance to the MSE definition in section V C With such a definition, residu 



als are expected to be zero mean and unitary variance in correspondence of the 'true' model Sa [fi] 0)- Different values 



of c were tested with a Monte Carlo estimation on the first 1000 data of the analysis reported in section V C Huber's 
norm is minimized for each MC iteration and the corresponding values of the parameter 9 are stored. The histogram 
of 6 shows a mean steadily fixed on 0.13 for values of c ranging from 0.001 to 0.1. Increasing the value of c, the 
distribution starts to shift toward the distribution obtained with MSE minimization. MSE and Huber distributions 
are practically undistinguishable for values c > 2. Since the 'true' value for the parameter 9 is set to 0.1, Huber's 
norm is performing better for values c < 0.1. It was then decided to use the value c = 0.05 for the final analysis 
reported in figure [Sj 
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